clear all
set more off
set mat 11000

set seed 913


version 10.0
mata:
mata set matalnum on


void coeffs1()
{
	real scalar rrrr


	st_view(V=.,.,tokens(st_local("base")))
	st_view(P=.,.,tokens(st_local("prone")))
	
	st_view(B=.,.,tokens(st_local("rowvar")))
	st_view(Z=.,.,tokens(st_local("instr")))
	st_view(N=.,.,st_local("cl"))
	rrrr=rows(B)*(cross(B,Z)/rows(B))*(cross(B,Z)/rows(B))'	
	st_matrix("A1", rrrr)
}


void coeffs2()
{
	numeric matrix II, g, g0, G 
	real scalar info

	st_view(S=.,.,tokens(st_local("score")))
	st_view(V=.,.,tokens(st_local("base")))
	st_view(P=.,.,tokens(st_local("prone")))
	
	st_view(B=.,.,tokens(st_local("rowvar")))
	st_view(Z=.,.,tokens(st_local("instr")))
	st_view(N=.,.,st_local("cl"))

	info=panelsetup(N,1)
	II=(Z:*B, V, P:*S, S)
	g=colsum(II)/rows(B)

	g0=J(cols(II),cols(II),0)	
	for (i=1; i<=rows(info); i++) {
		G=panelsubmatrix(II, i, info)			
		g0=colsum(G)'*colsum(G)+g0
	}

	g0=g0/rows(II)
	st_matrix("g1", g)
	st_matrix("A2", g0)
}



void coeffs3()
{
	numeric matrix II, g
	real scalar info

	st_view(S=.,.,tokens(st_local("score")))
	st_view(V=.,.,tokens(st_local("base")))
	st_view(P=.,.,tokens(st_local("prone")))
	
	st_view(B=.,.,tokens(st_local("rowvar")))
	st_view(Z=.,.,tokens(st_local("instr")))
	st_view(N=.,.,st_local("cl"))

	info=panelsetup(N,1)
	II=(Z:*B, V, P:*S, S)
	g=colsum(II)/rows(B)
	st_matrix("g1", g)
}

end

capture program quan drop
program quan, eclass
	version 10.0
	syntax varlist(max=1) [, end1(varname) end2(varname) instr(varlist) exog(varname) prone(varlist) tau(integer 50) ] 
	tempname temp cutoff below prb temp2 prb2 temp3 tempa belowa cutoffa temp4 temp5 temp4a temp5a temp6 resid moment ff temp4_t tmp temp3_t ff2
	local tau2=`tau'
	local tau=`tau'/100
	local cn=_N
	gen `ff'=1 if _n==1 
	egen `ff2'=sum(`ff')
	qui sum `ff2' 

	local cn=r(mean)
	local low1=9999999
	local high1=-9999999	
	local low2=9999999
	local high2=-9999999	


		
	tempname pdf cdf moment1 moment2 www temp5 temp3 temp5_t temp3_t moment2_t chi

	foreach var of varlist `prone' {
		tempname `var'_t
	}
	
	****initial estimate****
	local fff=1
	local num=-2000
	while `num'<=3500 {
	local num2=-2000
	while `num2'<=3500 {
		gen `temp'=`varlist'-`num'*`end1'-`num2'*`end2'
		egen `cutoff'=pctile(`temp'),  p(`tau2')			
		gen `below'=(`temp' <= `cutoff')
		
		qui logit `below' `prone'
		predict pr, rules
		gen `prb2'=pr	
		drop pr
		gen `temp5'=`prb2'-`below'
		local rowvar ""
		local rowvar  `"`rowvar' `temp5'"'
		local base = "`temp5'"	

	
		qui mata: coeffs1()
	
		local rrr2=A1[1,1]
		if (`fff'==1) {
			scalar residual1=`rrr2'				
			local beta1a=`num'
			local beta2a=`num2'
			local fff=2
		}
		else {

			if (`rrr2'<residual1) {
				local beta1a=`num'
				local beta2a=`num2'
				scalar residual1=`rrr2'

			}
		}
		drop `temp'-`temp5'
		local num2=`num2'+1
	}
		local num=`num'+1
	}


	***get variance matrix
	local num=`beta1a'
	local num2=`beta2a'

		gen `temp'=`varlist'-`num'*`end1'-`num2'*`end2'
				
		egen `cutoff'=pctile(`temp'),  p(`tau2')
		gen `below'=(`temp' <= `cutoff')

		qui logit `below' `prone'
		indeplist, local
		predict pr2, rules
		predict xb, xb rules
		gen `prb'=pr2
		drop pr2 
		
		gen `temp4'=`prb'-`below'

		local rowvar ""
		local rowvar  `"`rowvar' `temp4'"'

		gen `temp5'=`tau2'/100-`below'
	
		gen `moment2'=`below'-(exp(xb)/(1+exp(xb))) 
	
		drop xb		
		local base = "`temp5'"		
		local score = "`moment2'"

		qui mata: coeffs2()		
		matrix Sigma9=A2
			
			
		matrix R=g1*syminv(Sigma9)*g1'*`cn'
		matlist R
		local RRR2=R[1,1]

		drop `temp'-`moment2'




	local low1=`beta1a'
	local high1=`beta1a'	
	local low2=`beta2a'
	local high2=`beta2a'		

	******* INFERENCE ********

	
	**INFERENCE - low1
	local num=-2500
	local aaa=`low1'
	local hit=0
	while `num'<`aaa' {
	local num2=-2500

	local fff=1
	qui while `num2'<=3500 {
	if (`hit'==0) {
		gen `temp'=`varlist'-`num'*`end1'-`num2'*`end2'
				
		egen `cutoff'=pctile(`temp'),  p(`tau2')
		gen `below'=(`temp' <= `cutoff')

		qui logit `below' `prone'
		predict pr2, rules
		predict xb, xb rules

		gen `prb'=pr2
		drop pr2 

		gen `temp4'=`prb'-`below'
		local rowvar ""
		local rowvar  `"`rowvar' `temp4'"'

		gen `temp5'=`tau2'/100-`below'
		gen `moment2'=`below'-(exp(xb)/(1+exp(xb))) 
		drop xb		
		local base = "`temp5'"		
		local score = "`moment2'"

		qui mata: coeffs3()		
			
		matrix R=g1*syminv(Sigma9)*g1'*`cn'
		local RRR=R[1,1]-`RRR2'

		drop `temp'-`moment2'

		if (`RRR'<3.841) {
			local hit=1
			if (`num'<`low1') {
				local low1=`num'
			}
			if (`num'>`high1') {
				local high1=`num'
			}
			if (`num2'<`low2') {
				local low2=`num2'
			}
			if (`num2'>`high2') {
				local high2=`num2'
			}
		}	
	}
		local num2=`num2'+1	
	}
	

	local hit=0
	**INFERENCE - high1
	local num=3500
	while `num'>`high1' {
	local num2=-2500
	local fff=1
	qui while `num2'<=3500 {
	if (`hit'==0) {
		gen `temp'=`varlist'-`num'*`end1'-`num2'*`end2'
				
		egen `cutoff'=pctile(`temp'),  p(`tau2')
		gen `below'=(`temp' <= `cutoff')

		qui logit `below' `prone'
		predict pr2, rules
		predict xb, xb rules

		gen `prb'=pr2
		drop pr2 

		gen `temp4'=`prb'-`below'
		local rowvar ""
		local rowvar  `"`rowvar' `temp4'"'

		gen `temp5'=`tau2'/100-`below'
		gen `moment2'=`below'-(exp(xb)/(1+exp(xb))) 
		drop xb		
		local base = "`temp5'"		
		local score = "`moment2'"
	
		qui mata: coeffs3()		
	
		matrix R=g1*syminv(Sigma9)*g1'*`cn'
		local RRR=R[1,1]-`RRR2'

		drop `temp'-`moment2'

		if (`RRR'<3.841) {
			local hit=1
			if (`num'<`low1') {
				local low1=`num'
			}
			if (`num'>`high1') {
				local high1=`num'
			}
			if (`num2'<`low2') {
				local low2=`num2'
			}
			if (`num2'>`high2') {
				local high2=`num2'
			}
		}	

	}
		local num2=`num2'+1
	}

	local hit=0

	**INFERENCE - low2
	local num2=-2500
	local aaa=`low2'
	while `num2'<`aaa' {
	local num=-2500
	local fff=1
	qui while `num'<=3500 {
	if (`hit'==0) {
		gen `temp'=`varlist'-`num'*`end1'-`num2'*`end2'
				
		egen `cutoff'=pctile(`temp'),  p(`tau2')
		gen `below'=(`temp' <= `cutoff')

		qui logit `below' `prone'
		predict pr2, rules
		predict xb, xb rules

		gen `prb'=pr2
		drop pr2 

		gen `temp4'=`prb'-`below'
		local rowvar ""
		local rowvar  `"`rowvar' `temp4'"'

		gen `temp5'=`tau2'/100-`below'
		gen `moment2'=`below'-(exp(xb)/(1+exp(xb))) 
		drop xb		
		local base = "`temp5'"		
		local score = "`moment2'"
	
		qui mata: coeffs3()		
		
				
		matrix R=g1*syminv(Sigma9)*g1'*`cn'
		local RRR=R[1,1]-`RRR2'

		drop `temp'-`moment2'

		if (`RRR'<3.841) {
			local hit=1
			if (`num'<`low1') {
				local low1=`num'
			}
			if (`num'>`high1') {
				local high1=`num'
			}
			if (`num2'<`low2') {
				local low2=`num2'
			}
			if (`num2'>`high2') {
				local high2=`num2'
			}
		}	

	}
		local num=`num'+1
	}



	local hit=0


	**INFERENCE - high2
	local num2=3500
	while `num2'>`high2' {
	local num=-2500
	local fff=1
	qui while `num'<=3500 {
	if (`hit'==0) {
		gen `temp'=`varlist'-`num'*`end1'-`num2'*`end2'
				
		egen `cutoff'=pctile(`temp'),  p(`tau2')
		gen `below'=(`temp' <= `cutoff')


		qui logit `below' `prone'
		predict pr2, rules
		predict xb, xb rules


		gen `prb'=pr2
		drop pr2 

		gen `temp4'=`prb'-`below'
		local rowvar ""
		local rowvar  `"`rowvar' `temp4'"'

		gen `temp5'=`tau2'/100-`below'
		gen `moment2'=`below'-(exp(xb)/(1+exp(xb))) 
		drop xb		
		local base = "`temp5'"		
		local score = "`moment2'"
	
		qui mata: coeffs3()		
		
				
		matrix R=g1*syminv(Sigma9)*g1'*`cn'
		local RRR=R[1,1]-`RRR2'

		drop `temp'-`moment2'

		if (`RRR'<3.841) {
			local hit=1
			if (`num'<`low1') {
				local low1=`num'
			}
			if (`num'>`high1') {
				local high1=`num'
			}
			if (`num2'<`low2') {
				local low2=`num2'
			}
			if (`num2'>`high2') {
				local high2=`num2'
			}
		}	
	}
		local num=`num'+1
	}


	
	ereturn scalar low1=`low1'
	ereturn scalar high1=`high1'
	ereturn scalar low2=`low2'
	ereturn scalar high2=`high2'
	ereturn scalar b2=`beta2a'
	ereturn scalar b1=`beta1a'

	
end
